###########################################
#
# APPENDIX for Habers and Ingram PSRM 2018
# Paper: "Spatial Tools for Case Selection: Using LISAs to Design Mixed-Methods Research"
#
###########################################
#
# NOTE:
# This file is designed to be run from main R script for replication.
# See README file.
# To run this file from main script, uncomment line near end of main 
# script that sources this file (e.g., "source()")
#
# To run this file on its own, apart from main script, then need to
# uncomment all lines below within next three sections, including:
# (1) INSTALL PACKAGES
# (2) SET WORKING DIRECTORY AND SUBDIRECTORIES
# (3) LOAD DATA
#
###########################################

###############################
# INSTALL PACKAGES
###############################

# install.packages("pacman")
#library(pacman)
#pacman::p_load(spdep, graphicsQC, prettyR, RColorBrewer, GWmodel, car, spgwr, grid, 
#               gridExtra, classInt, rgeos, maptools, 
#               RgoogleMaps, rgdal, weights, xts, spacetime, plm, 
#               ggplot2, maps, plm, plyr, 
#               parallel, scales, cowplot, AER, 
#               foreign, haven, readstata13, psych)

###############################################
# SET WORKING DIRECTORY AND SUBDIRECTORIES
###############################################

# note: this is set to my own working directory; select your own

#path <- 'C:/Users/mi122167/Dropbox/Met Matt/PSRM2018/replicationmaterials2018'
#setwd(path)

# create subdirectories
#dir.create('./code', showWarnings = TRUE)
#dir.create('./data', showWarnings = TRUE)
#dir.create('./data/original', showWarnings = TRUE)
#dir.create('./data/working', showWarnings = TRUE)
#dir.create('./figures', showWarnings = TRUE)
#dir.create('./shapefiles', showWarnings = TRUE)

##############################
# LOAD DATA
##############################
# if returning to data after analysis, just load workspace (.RData file)
# if no data yet, then skip this line
#load("./data/working/harbersingram_psrm_20181109.RData")


#########################################################
# APPENDIX CONTENT
#########################################################
#
# color version of Figure 1; reported in appendix
#
#########################################################

DV <- shp$HR90
quadrant <- vector(mode="numeric",length=nrow(shp@data))
# center on mean
cDV <- DV - mean(DV) 
lagDV <- lag.listw(w, DV)
# center on mean
clagDV <- lagDV - mean(lagDV)

cl <- as.factor(shp$ls.py.clper)
sig <- shp$ls.py.pperm

# fig1a - moran plot y=DV, color
dat2 = data.frame(y = cDV, Wy = clagDV, cl=cl, sig=sig, meany=mean(cDV), meanWy=mean(clagDV))

g <- ggplot(dat2, aes(x = y, y = Wy))+
  geom_point(colour="black", pch=21, size=3,
             aes(fill = factor(cl))) +
  scale_fill_manual(name = "Cluster",
                    values = c("0" = "white",
                               "1" = "red",
                               "2" = "blue",
                               "3" = "lightblue",
                               "4" = "pink"),
                    labels = c("n.s.", "high-high", "low-low", "low-high", "high-low")) +
  #labs(x="", y="", title="LISA Clusters (DV: HR90)") +
  geom_smooth(method = "lm", se=F, colour="black", size=.7) + 
  geom_vline(xintercept=dat2$meany,colour="black",linetype="longdash") + 
  geom_hline(yintercept=dat2$meanWy,colour="black",linetype="longdash")+ 
  xlab("HR90 (centered on mean)") + ylab("spatial lag HR90 (centered on mean)")+
  theme(axis.line=element_line(color="black"),
        axis.title.x=element_text(size=10,vjust=0.1),
        axis.title.y=element_text(size=10,vjust=0.1),
        axis.text= element_text(colour="black", size=10, angle=0,face = "plain"),
        #plot.title = element_text(size = 10, lineheight=.8, face="bold", vjust=1), # make title bold and add space
        legend.title = element_text(size = 8), # legend title size
        legend.text = element_text(size = 8), # legend label size
        panel.grid.major = element_blank(), panel.grid.minor = element_blank(), # removes background grid
        panel.background = element_blank(), # removes grey background; could add black axis lines with #axis.line = element_line(colour = "black")) 
        panel.spacing=unit(c(0,0,0,0), "lines"),
        plot.margin=unit(c(0,0,0,0), "mm"))  # sets margin around full plot at top, right, bottom, and left; units can also be "lines" or "cm"

tiff(file="./figures/appendix/fig1A_moranploty_col_cen.tif", height=6, width=6, units="in", res=300) # could also do pdf, jpeg, bmp, tiff
print(g)
dev.off()


#################################################
#
# Color version of Figure 2
#
#################################################

# load shapefile with all data from python analysis (as well as GeoDa)
shp2 <- readOGR(dsn='./shapefiles', layer="NAT3-8")

shp <- shp2

shp@data$id = rownames(shp@data)
shp.points <- fortify(shp, region="id")

# if get error about gpclibPermit not being TRUE, try re-installing rgeos:
#install.packages("rgeos")
# if that does not work, try:
#install.packages("gpclib", type="source")
#library(gpclib)
#gpclibPermit()

shp.df = join(shp.points, shp@data, by="id")
#shp.df[shp.df==0] <- NA

#names(shp.df)

# Figure 2 raw HR90 values

plot1.1 <- ggplot(shp.df) + 
  aes(long,lat,group=group,fill=HR90) + 
  geom_polygon(colour="transparent", fill="white")
#plot1.1

plot1.1a <- plot1.1 + 
  geom_polygon(aes(x = long, y = lat, group = group), data = shp.df) +            
  #             colour = 'white', fill = 'black', alpha = .4, size = .3) +
  #scale_fill_distiller(name="LISA values", palette = "Reds", 
  #                    trans="reverse", 
  # breaks = pretty_breaks(n = 5), 
  #                   na.value="white", guide="colourbar") + 
  #scale_fill_gradient(name="LISA values", values=c("white", "red", "blue","lightblue","pink"), breaks = c("0", "1", "2", "3", "4"), labels=c("n.s.", "high-high", "low-low","low-high","hgh-low")) + 
  scale_fill_gradientn(name="HR90", 
                       #values=rescale(c(0, 20, 60)), 
                       colours=c("white","red")) +
  labs(x="", y="", title="Homicide Rates (HR90)") +
  theme(axis.ticks.y = element_blank(),axis.text.y = element_blank(), # get rid of x ticks/text
        axis.ticks.x = element_blank(),axis.text.x = element_blank(), # get rid of y ticks/text
        axis.line.x = element_blank(), axis.line.y = element_blank(), # get rid of axis lines, if present
        plot.title = element_text(size = 10, lineheight=.8, face="bold", vjust=1), # make title bold and add space
        legend.title = element_text(size = 8), # legend title size
        legend.text = element_text(size = 8), # legend label size
        panel.grid.major = element_blank(), panel.grid.minor = element_blank(), # removes background grid
        panel.background = element_blank(), # removes grey background; could add black axis lines with #axis.line = element_line(colour = "black")) 
        panel.spacing=unit(c(0,0,0,0), "lines"),
        plot.margin=unit(c(0,0,0,0), "mm")) + # sets margin around full plot at top, right, bottom, and left; units can also be "lines" or "cm"
  coord_equal(ratio=1) +
  geom_polygon(data=shp.st, aes(long,lat, group=group), colour="grey50", fill=NA) 
#plot1.1a

#tiff(file="./figures/figA_hr90values-color.tif", height=6, width=9, units="in", res=300) # could also do pdf, jpeg, bmp, tiff
#print(plot1.1a)
#dev.off()

#####################################
# Figure 2 - color
#####################################

# LISA values map - color

plot1p <- ggplot(shp.df) + 
  aes(long,lat,group=group,fill=ls.py.I) + 
  geom_polygon(colour="transparent", fill="white")
#plot1p

plot1pa <- plot1p + 
  geom_polygon(aes(x = long, y = lat, group = group), data = shp.df) +            
  #             colour = 'white', fill = 'black', alpha = .4, size = .3) +
  #scale_fill_distiller(name="LISA values", palette = "Reds", 
  #                    trans="reverse", 
  # breaks = pretty_breaks(n = 5), 
  #                   na.value="white", guide="colourbar") + 
  #scale_fill_gradient(name="LISA values", values=c("white", "red", "blue","lightblue","pink"), breaks = c("0", "1", "2", "3", "4"), labels=c("n.s.", "high-high", "low-low","low-high","hgh-low")) + 
  #scale_fill_gradientn(name="LISA values", values=rescale(c(12, 4, 0, -2, -6)), colours=c("red","white","blue")) +
  scale_fill_gradient2(name="LISA values", low="blue", mid="white", high="red", midpoint=0) +
  labs(x="", y="", title="LISA Values (DV: HR90)") +
  theme(axis.ticks.y = element_blank(),axis.text.y = element_blank(), # get rid of x ticks/text
        axis.ticks.x = element_blank(),axis.text.x = element_blank(), # get rid of y ticks/text
        axis.line.x = element_blank(), axis.line.y = element_blank(), # get rid of axis lines, if present
        plot.title = element_text(size = 10, lineheight=.8, face="bold", vjust=1), # make title bold and add space
        legend.title = element_text(size = 8), # legend title size
        legend.text = element_text(size = 8), # legend label size
        panel.grid.major = element_blank(), panel.grid.minor = element_blank(), # removes background grid
        panel.background = element_blank(), # removes grey background; could add black axis lines with #axis.line = element_line(colour = "black")) 
        panel.spacing=unit(c(0,0,0,0), "lines"),
        plot.margin=unit(c(0,0,0,0), "mm")) + # sets margin around full plot at top, right, bottom, and left; units can also be "lines" or "cm"
  coord_equal(ratio=1) +
  geom_polygon(data=shp.st, aes(long,lat, group=group), colour="grey50", fill=NA) 
#plot1pa

#tiff(file="./figures/figA_lisavalues-perm-color.tif", height=6, width=9, units="in", res=300) # could also do pdf, jpeg, bmp, tiff
#print(plot1pa)
#dev.off()

# Map of LISA p-values - color

plot2p <- ggplot(shp.df) + 
  aes(long,lat,group=group,fill=ls.py.pperm) + 
  geom_polygon(colour="transparent", fill="white")
#plot2p

plot2pa <- plot2p + 
  geom_polygon(aes(x = long, y = lat, group = group), data = subset(shp.df, ls.py.pperm<.2)) +            
  #             colour = 'white', fill = 'black', alpha = .4, size = .3) +
  scale_fill_gradientn(name="LISA p", 
                       #values=rescale(c(0, .05, .2)), 
                       colours=c("darkred", "white"), guide="colourbar") + 
  #scale_fill_gradient(name="LISA sig (p)", low="white", high="darkred", guide="colourbar") + 
  labs(x="", y="", title="LISA p-values (DV: HR90)") +
  theme(axis.ticks.y = element_blank(),axis.text.y = element_blank(), # get rid of x ticks/text
        axis.ticks.x = element_blank(),axis.text.x = element_blank(), # get rid of y ticks/text
        axis.line.x = element_blank(), axis.line.y = element_blank(), # get rid of axis lines, if present
        plot.title = element_text(size = 10, lineheight=.8, face="bold", vjust=1), # make title bold and add space
        legend.title = element_text(size = 8), # legend title size
        legend.text = element_text(size = 8), # legend label size
        panel.grid.major = element_blank(), panel.grid.minor = element_blank(), # removes background grid
        panel.background = element_blank(), # removes grey background; could add black axis lines with #axis.line = element_line(colour = "black")) 
        panel.spacing=unit(c(0,0,0,0), "lines"),
        plot.margin=unit(c(0,0,0,0), "mm")) + # sets margin around full plot at top, right, bottom, and left; units can also be "lines" or "cm"
  coord_equal(ratio=1) +
  geom_polygon(data=shp.st, aes(long,lat, group=group), colour="grey50", fill=NA) 
#plot2pa

#tiff(file="./figures/figA_lisap-perm.tif", height=6, width=9, units="in", res=300) # could also do pdf, jpeg, bmp, tiff
#print(plot2pa)
#dev.off()

# Map of LISA clusters with permutation approach in python - color

plot3pyperm <- ggplot(shp.df) + 
  aes(long,lat,group=group,fill=as.factor(ls.py.clper)) + 
  geom_polygon(colour="transparent", fill="white")
#plot3perm

plot3pyperma <- plot3pyperm + 
  geom_polygon(aes(x = long, y = lat, group = group), data = shp.df) +            
  #             colour = 'white', fill = 'black', alpha = .4, size = .3) +
  #scale_fill_distiller(name=NULL, palette = "Greys", trans="reverse", breaks = pretty_breaks(n = 5), na.value="white", guide="colourbar") + 
  scale_fill_manual(name="LISA cluster", values=c("white", "red", "blue","lightblue","pink"), breaks = c("0", "1", "2", "3", "4"), 
                    labels=c("n.s.", "high-high", "low-low","low-high","high-low"), guide="legend") + 
  labs(x="", y="", title="LISA Clusters (DV: HR90)") +
  theme(axis.ticks.y = element_blank(),axis.text.y = element_blank(), # get rid of x ticks/text
        axis.ticks.x = element_blank(),axis.text.x = element_blank(), # get rid of y ticks/text
        axis.line.x = element_blank(), axis.line.y = element_blank(), # get rid of axis lines, if present
        plot.title = element_text(size = 10, lineheight=.8, face="bold", vjust=1), # make title bold and add space
        legend.title = element_text(size = 8), # legend title size
        legend.text = element_text(size = 8), # legend label size
        panel.grid.major = element_blank(), panel.grid.minor = element_blank(), # removes background grid
        panel.background = element_blank(), # removes grey background; could add black axis lines with #axis.line = element_line(colour = "black")) 
        panel.spacing=unit(c(0,0,0,0), "lines"),
        plot.margin=unit(c(0,0,0,0), "mm")) + # sets margin around full plot at top, right, bottom, and left; units can also be "lines" or "cm"
  coord_equal(ratio=1) +
  geom_polygon(data=shp.st, aes(long,lat, group=group), colour="grey50", fill=NA) 
# 
#plot3pyperma

#tiff(file="./figures/figA_lisa-permpy-color.tif", height=6, width=9, units="in", res=300) # could also do pdf, jpeg, bmp, tiff
#print(plot3pyperma)
#dev.off()

###################################################
# Combined Graphs with Permutation Approach - COLOR
###################################################

# DV (Figure 2)

tiff(file="./figures/appendix/fig2A_comb_pyperm_color.tif", height=9, width=6, units="in", res=300) # could also do pdf, jpeg, bmp, tiff
grid.arrange(plot1.1a, plot1pa, plot2pa, plot3pyperma, ncol=1)
dev.off()


##########################################################
#
# Alternative saddlepoint estimation version of Figure 2
#
##########################################################

nb <- w_10nnb
lisa_s <- as.data.frame(summary(localmoran.sad(lm(HR90 ~ 1, shp), nb=nb, style="C"))) # this works; need summary for saddlepoint syntax

# create LISA cluster identifiers
DV <- shp$HR90
quadrant <- vector(mode="numeric",length=nrow(lisa_s))
cDV <- DV - mean(DV) 
lagDV <- lag.listw(w, DV)
clagDV <- lagDV - mean(lagDV)

# LISA significance with saddlepoint method
p <- lisa_s[,5]
quadrant <- vector(mode="numeric",length=nrow(lisa_s))
quadrant[cDV >0 & clagDV>0 & p<=.05] <- 1 
quadrant[cDV <0 & clagDV<0 & p<=.05] <- 2      
quadrant[cDV <0 & clagDV>0 & p<=.05] <- 3
quadrant[cDV >0 & clagDV<0 & p<=.05] <- 4
# non-significant obs will remain coded as zeroes (0)

# merge LISA of DV data into shapefile
# merge a single variable from table
shp$lisa_s.dv <- lisa_s[,4]
shp$lisa_s.p.dv <- p
shp$lisa_s.cl.dv <- as.factor(quadrant)
#names(shp)


shp@data$id = rownames(shp@data)
shp.points <- fortify(shp, region="id")

# if get error about gpclibPermit not being TRUE, try re-installing rgeos:
#install.packages("rgeos")
# if that does not work, try:
#install.packages("gpclib", type="source")
#library(gpclib)
#gpclibPermit()

shp.df = join(shp.points, shp@data, by="id")
#shp.df[shp.df==0] <- NA


# Figures
plot1s <- ggplot(shp.df) + 
  aes(long,lat,group=group,fill=lisa_s.dv) + 
  geom_polygon(colour="transparent", fill="white")
#plot1

plot1sa <- plot1s + 
  geom_polygon(aes(x = long, y = lat, group = group), data = shp.df) +            
  #             colour = 'white', fill = 'black', alpha = .4, size = .3) +
  #scale_fill_distiller(name="LISA values", palette = "Reds", 
  #                    trans="reverse", 
  # breaks = pretty_breaks(n = 5), 
  #                   na.value="white", guide="colourbar") + 
  #scale_fill_gradient(name="LISA values", values=c("white", "red", "blue","lightblue","pink"), breaks = c("0", "1", "2", "3", "4"), labels=c("n.s.", "high-high", "low-low","low-high","hgh-low")) + 
  scale_fill_gradientn(name="LISA values", 
                       #values=rescale(c(12, 4, 0, -2, -6)), 
                       colours=c("red","white","blue")) +
  labs(x="", y="", title="LISA Values (DV: HR90)") +
  theme(axis.ticks.y = element_blank(),axis.text.y = element_blank(), # get rid of x ticks/text
        axis.ticks.x = element_blank(),axis.text.x = element_blank(), # get rid of y ticks/text
        axis.line.x = element_blank(), axis.line.y = element_blank(), # get rid of axis lines, if present
        plot.title = element_text(size = 10, lineheight=.8, face="bold", vjust=1), # make title bold and add space
        legend.title = element_text(size = 8), # legend title size
        legend.text = element_text(size = 8), # legend label size
        panel.grid.major = element_blank(), panel.grid.minor = element_blank(), # removes background grid
        panel.background = element_blank(), # removes grey background; could add black axis lines with #axis.line = element_line(colour = "black")) 
        panel.spacing=unit(c(0,0,0,0), "lines"),
        plot.margin=unit(c(0,0,0,0), "mm")) + # sets margin around full plot at top, right, bottom, and left; units can also be "lines" or "cm"
  coord_equal(ratio=1) +
  geom_polygon(data=shp.st, aes(long,lat, group=group), colour="grey50", fill=NA) 
#plot1sa

#tiff(file="./figures/figA_lisavalues_sad_color.tif", height=6, width=9, units="in", res=300) # could also do pdf, jpeg, bmp, tiff
#print(plot1sa)
#dev.off()


# LISA p-values with saddlepoint method

plot2s <- ggplot(shp.df) + 
  aes(long,lat,group=group,fill=lisa_s.p.dv) + 
  geom_polygon(colour="transparent", fill="white")
#plot2s

plot2sa <- plot2s + 
  geom_polygon(aes(x = long, y = lat, group = group), data = subset(shp.df, lisa_s.p.dv<.2)) +            
  #             colour = 'white', fill = 'black', alpha = .4, size = .3) +
  scale_fill_gradientn(name="LISA p", 
                       #values=rescale(c(0, .05, .2)), 
                       colours=c("darkred", "white"), guide="colourbar") + 
  #scale_fill_gradient(name="LISA sig (p)", low="white", high="darkred", guide="colourbar") + 
  labs(x="", y="", title="LISA p-values (DV: HR90)") +
  theme(axis.ticks.y = element_blank(),axis.text.y = element_blank(), # get rid of x ticks/text
        axis.ticks.x = element_blank(),axis.text.x = element_blank(), # get rid of y ticks/text
        axis.line.x = element_blank(), axis.line.y = element_blank(), # get rid of axis lines, if present
        plot.title = element_text(size = 10, lineheight=.8, face="bold", vjust=1), # make title bold and add space
        legend.title = element_text(size = 8), # legend title size
        legend.text = element_text(size = 8), # legend label size
        panel.grid.major = element_blank(), panel.grid.minor = element_blank(), # removes background grid
        panel.background = element_blank(), # removes grey background; could add black axis lines with #axis.line = element_line(colour = "black")) 
        panel.spacing=unit(c(0,0,0,0), "lines"),
        plot.margin=unit(c(0,0,0,0), "mm")) + # sets margin around full plot at top, right, bottom, and left; units can also be "lines" or "cm"
  coord_equal(ratio=1) +
  geom_polygon(data=shp.st, aes(long,lat, group=group), colour="grey50", fill=NA) 
#plot2sa

#tiff(file="./figures/figA_lisap_sad-col.tif", height=6, width=9, units="in", res=300) # could also do pdf, jpeg, bmp, tiff
#print(plot2sa)
#dev.off()

# black and white p-value with saddlepoint

plot2sb <- plot2s + 
  geom_polygon(aes(x = long, y = lat, group = group), data = subset(shp.df, lisa_s.p.dv<.2)) +            
  #             colour = 'white', fill = 'black', alpha = .4, size = .3) +
  scale_fill_gradientn(name="LISA p", 
                       #values=rescale(c(0, .05, .2)), 
                       colours=c("black", "white"), guide="colourbar") + 
  #scale_fill_gradient(name="LISA sig (p)", low="white", high="darkred", guide="colourbar") + 
  labs(x="", y="", title="LISA p-values (DV: HR90)") +
  theme(axis.ticks.y = element_blank(),axis.text.y = element_blank(), # get rid of x ticks/text
        axis.ticks.x = element_blank(),axis.text.x = element_blank(), # get rid of y ticks/text
        axis.line.x = element_blank(), axis.line.y = element_blank(), # get rid of axis lines, if present
        plot.title = element_text(size = 10, lineheight=.8, face="bold", vjust=1), # make title bold and add space
        legend.title = element_text(size = 8), # legend title size
        legend.text = element_text(size = 8), # legend label size
        panel.grid.major = element_blank(), panel.grid.minor = element_blank(), # removes background grid
        panel.background = element_blank(), # removes grey background; could add black axis lines with #axis.line = element_line(colour = "black")) 
        panel.spacing=unit(c(0,0,0,0), "lines"),
        plot.margin=unit(c(0,0,0,0), "mm")) + # sets margin around full plot at top, right, bottom, and left; units can also be "lines" or "cm"
  coord_equal(ratio=1) +
  geom_polygon(data=shp.st, aes(long,lat, group=group), colour="grey50", fill=NA) 
#plot2sb

#tiff(file="./figures/figA_lisap_sad_bw.tif", height=6, width=9, units="in", res=300) # could also do pdf, jpeg, bmp, tiff
#print(plot2sb)
#dev.off()


# LISA cluster with saddlepoint method

plot3s <- ggplot(shp.df) + 
  aes(long,lat,group=group,fill=lisa_s.cl.dv) + 
  geom_polygon(colour="transparent", fill="white")
#plot3s

plot3sa <- plot3s + 
  geom_polygon(aes(x = long, y = lat, group = group), data = shp.df) +            
  #             colour = 'white', fill = 'black', alpha = .4, size = .3) +
  #scale_fill_distiller(name=NULL, palette = "Greys", trans="reverse", breaks = pretty_breaks(n = 5), na.value="white", guide="colourbar") + 
  scale_fill_manual(name="LISA cluster", 
                    values=c("white", "red", "blue","lightblue","pink"), 
                    breaks = c("0", "1", "2", "3", "4"), 
                    labels=c("n.s.", "high-high", "low-low","low-high","hgh-low"), guide="legend") + 
  labs(x="", y="", title="LISA Clusters (DV: HR90)") +
  theme(axis.ticks.y = element_blank(),axis.text.y = element_blank(), # get rid of x ticks/text
        axis.ticks.x = element_blank(),axis.text.x = element_blank(), # get rid of y ticks/text
        axis.line.x = element_blank(), axis.line.y = element_blank(), # get rid of axis lines, if present
        plot.title = element_text(size = 10, lineheight=.8, face="bold", vjust=1), # make title bold and add space
        legend.title = element_text(size = 8), # legend title size
        legend.text = element_text(size = 8), # legend label size
        panel.grid.major = element_blank(), panel.grid.minor = element_blank(), # removes background grid
        panel.background = element_blank(), # removes grey background; could add black axis lines with #axis.line = element_line(colour = "black")) 
        panel.spacing=unit(c(0,0,0,0), "lines"),
        plot.margin=unit(c(0,0,0,0), "mm")) + # sets margin around full plot at top, right, bottom, and left; units can also be "lines" or "cm"
  coord_equal(ratio=1) +
  geom_polygon(data=shp.st, aes(long,lat, group=group), colour="grey50", fill=NA) 
#plot3sa

#tiff(file="./figures/figA_lisa-sad-color.tif", height=6, width=9, units="in", res=300) # could also do pdf, jpeg, bmp, tiff
#print(plot3sa)
#dev.off()

# black and white clusters with saddlepoint

plot3sb <- plot3s + 
  geom_polygon(aes(x = long, y = lat, group = group), data = shp.df) +            
  #             colour = 'white', fill = 'black', alpha = .4, size = .3) +
  #scale_fill_distiller(name=NULL, palette = "Greys", trans="reverse", breaks = pretty_breaks(n = 5), na.value="white", guide="colourbar") + 
  scale_fill_manual(name="LISA cluster", values=c("white", "black", "grey","",""), breaks = c("0", "1", "2", "3", "4"), labels=c("n.s.", "high-high", "low-low","low-high","hgh-low")) + 
  labs(x="", y="", title="LISA Clusters (DV: HR90)") +
  theme(axis.ticks.y = element_blank(),axis.text.y = element_blank(), # get rid of x ticks/text
        axis.ticks.x = element_blank(),axis.text.x = element_blank(), # get rid of y ticks/text
        axis.line.x = element_blank(), axis.line.y = element_blank(), # get rid of axis lines, if present
        plot.title = element_text(size = 10, lineheight=.8, face="bold", vjust=1), # make title bold and add space
        legend.title = element_text(size = 8), # legend title size
        legend.text = element_text(size = 8), # legend label size
        panel.grid.major = element_blank(), panel.grid.minor = element_blank(), # removes background grid
        panel.background = element_blank(), # removes grey background; could add black axis lines with #axis.line = element_line(colour = "black")) 
        panel.spacing=unit(c(0,0,0,0), "lines"),
        plot.margin=unit(c(0,0,0,0), "mm")) + # sets margin around full plot at top, right, bottom, and left; units can also be "lines" or "cm"
  coord_equal(ratio=1) +
  geom_polygon(data=shp.st, aes(long,lat, group=group), colour="grey50", fill=NA) 
#plot3sb

#tiff(file="./figures/figA_lisa-bw.tif", height=6, width=9, units="in", res=300) # could also do pdf, jpeg, bmp, tiff
#print(plot3sb)
#dev.off()


#################################################################
#
# COMBINED GRAPH OF FIGURE 2 (DV) with SADDLEPOINT ESTIMATION
#
#################################################################

tiff(file="./figures/appendix/fig2A_comb_dv-sad_color.tif", height=9, width=6, units="in", res=300) # could also do pdf, jpeg, bmp, tiff
grid.arrange(plot1.1a, plot1sa, plot2sa, plot3sa, ncol=1)
dev.off()


####################################################
### Moran plot of residuals with ggplot2 - color
####################################################

res <- ols2$residuals
# center on mean
lagres <- lag.listw(w, res)
# center on mean
#clagDV <- lagDV - mean(lagDV)

cl <- as.factor(shp$ls.py.qresp)

dat2 = data.frame(res = res, Wres = lagres, cl=cl, meanres=mean(res), meanWres=mean(lagres))


g <- ggplot(dat2, aes(x = res, y = Wres))+
  geom_point(colour="black", pch=21, size=3,
             aes(fill = factor(cl))) +
  scale_fill_manual(name = "Cluster",
                    values = c("0" = "white",
                               "1" = "red",
                               "2" = "blue",
                               "3" = "lightblue",
                               "4" = "pink"),
                    labels = c("n.s.", "high-high", "low-low", "low-high", "high-low")) +
  #labs(x="", y="", title="LISA Clusters (DV: HR90)") +
  geom_smooth(method = "lm", se=F, colour="black", size=.7) + 
  geom_vline(xintercept=dat2$meanres,colour="black",linetype="longdash") + 
  geom_hline(yintercept=dat2$meanWres,colour="black",linetype="longdash")+ 
  xlab("residual (centered on mean)") + ylab("spatial lag residual (centered on mean)")+
  theme(axis.line=element_line(color="black"),
        axis.title.x=element_text(size=10,vjust=0.1),
        axis.title.y=element_text(size=10,vjust=0.1),
        axis.text= element_text(colour="black", size=10, angle=0,face = "plain"),
        #plot.title = element_text(size = 10, lineheight=.8, face="bold", vjust=1), # make title bold and add space
        legend.title = element_text(size = 8), # legend title size
        legend.text = element_text(size = 8), # legend label size
        panel.grid.major = element_blank(), panel.grid.minor = element_blank(), # removes background grid
        panel.background = element_blank(), # removes grey background; could add black axis lines with #axis.line = element_line(colour = "black")) 
        panel.spacing=unit(c(0,0,0,0), "lines"),
        plot.margin=unit(c(0,0,0,0), "mm"))  # sets margin around full plot at top, right, bottom, and left; units can also be "lines" or "cm"

tiff(file="./figures/appendix/fig1A_moranplotres_col.tif", height=6, width=6, units="in", res=300) # could also do pdf, jpeg, bmp, tiff
print(g)
dev.off()



########################################3
#
# LISA maps of residuals in color
#
#########################################

# Residuals (ols2res)

plot4 <- ggplot(shp.df) + 
  aes(long,lat,group=group,fill=ols2res) + 
  geom_polygon(colour="transparent", fill="white")
#plot4

plot4a <- plot4 + 
  geom_polygon(aes(x = long, y = lat, group = group), data = shp.df) +            
  #             colour = 'white', fill = 'black', alpha = .4, size = .3) +
  #scale_fill_distiller(name="LISA values", palette = "Reds", 
  #                    trans="reverse", 
  # breaks = pretty_breaks(n = 5), 
  #                   na.value="white", guide="colourbar") + 
  #scale_fill_gradient(name="LISA values", values=c("white", "red", "blue","lightblue","pink"), breaks = c("0", "1", "2", "3", "4"), labels=c("n.s.", "high-high", "low-low","low-high","hgh-low")) + 
  #scale_fill_gradientn(name="Residuals", values=rescale(c(69, 4, 0, -2, -20)), colours=c("red","white","blue")) +
  scale_fill_gradient2(name="Residuals", high="red",mid="white",low="blue", midpoint=0) +
  labs(x="", y="", title="OLS Residuals)") +
  theme(axis.ticks.y = element_blank(),axis.text.y = element_blank(), # get rid of x ticks/text
        axis.ticks.x = element_blank(),axis.text.x = element_blank(), # get rid of y ticks/text
        axis.line.x = element_blank(), axis.line.y = element_blank(), # get rid of axis lines, if present
        plot.title = element_text(size = 10, lineheight=.8, face="bold", vjust=1), # make title bold and add space
        legend.title = element_text(size = 8), # legend title size
        legend.text = element_text(size = 8), # legend label size
        panel.grid.major = element_blank(), panel.grid.minor = element_blank(), # removes background grid
        panel.background = element_blank(), # removes grey background; could add black axis lines with #axis.line = element_line(colour = "black")) 
        panel.spacing=unit(c(0,0,0,0), "lines"),
        plot.margin=unit(c(0,0,0,0), "mm")) + # sets margin around full plot at top, right, bottom, and left; units can also be "lines" or "cm"
  coord_equal(ratio=1) +
  geom_polygon(data=shp.st, aes(long,lat, group=group), colour="grey50", fill=NA) 
#plot4a

#tiff(file="./figures/figA_olsresiduals-col.tif", height=6, width=9, units="in", res=300) # could also do pdf, jpeg, bmp, tiff
#print(plot4a)
#dev.off()

# lisa values of residuals using permutation approach in python - color

plot4.2pyperm <- ggplot(shp.df) + 
  aes(long,lat,group=group,fill=ls.py.lresp) + 
  geom_polygon(colour="transparent", fill="white")
#plot4.2pyperm

plot4bpyperma <- plot4.2pyperm + 
  geom_polygon(aes(x = long, y = lat, group = group), data = shp.df) +            
  #             colour = 'white', fill = 'black', alpha = .4, size = .3) +
  #scale_fill_distiller(name="LISA values", palette = "Reds", 
  #                    trans="reverse", 
  # breaks = pretty_breaks(n = 5), 
  #                   na.value="white", guide="colourbar") + 
  #scale_fill_gradient(name="LISA values", values=c("white", "red", "blue","lightblue","pink"), breaks = c("0", "1", "2", "3", "4"), labels=c("n.s.", "high-high", "low-low","low-high","hgh-low")) + 
  #scale_fill_gradientn(name="LISA values", values=rescale(c(7, 0, -7)), colours=c("red","white","blue")) +
  scale_fill_gradient2(name="LISA values", high="red", mid="white", low="blue", midpoint=0) +
  labs(x="", y="", title="LISA Values (residuals)") +
  theme(axis.ticks.y = element_blank(),axis.text.y = element_blank(), # get rid of x ticks/text
        axis.ticks.x = element_blank(),axis.text.x = element_blank(), # get rid of y ticks/text
        axis.line.x = element_blank(), axis.line.y = element_blank(), # get rid of axis lines, if present
        plot.title = element_text(size = 10, lineheight=.8, face="bold", vjust=1), # make title bold and add space
        legend.title = element_text(size = 8), # legend title size
        legend.text = element_text(size = 8), # legend label size
        panel.grid.major = element_blank(), panel.grid.minor = element_blank(), # removes background grid
        panel.background = element_blank(), # removes grey background; could add black axis lines with #axis.line = element_line(colour = "black")) 
        panel.spacing=unit(c(0,0,0,0), "lines"),
        plot.margin=unit(c(0,0,0,0), "mm")) + # sets margin around full plot at top, right, bottom, and left; units can also be "lines" or "cm"
  coord_equal(ratio=1) +
  geom_polygon(data=shp.st, aes(long,lat, group=group), colour="grey50", fill=NA) 
plot4bpyperma


#tiff(file="./figures/figA_lisavalues-pyperm-res-color.tif", height=6, width=9, units="in", res=300) # could also do pdf, jpeg, bmp, tiff
#print(plot4bpyperma)
#dev.off()

# LISA p-values permutation approach in python - color

plot5pyperm <- ggplot(shp.df) + 
  aes(long,lat,group=group,fill=ls.py.presp) + 
  geom_polygon(colour="transparent", fill="white")
#plot5pyperm

plot5pyperma <- plot5pyperm + 
  geom_polygon(aes(x = long, y = lat, group = group), data = subset(shp.df, ls.py.presp<.2)) +            
  #             colour = 'white', fill = 'black', alpha = .4, size = .3) +
  scale_fill_gradientn(name="LISA p", 
                       #values=rescale(c(0, .05, .2)), 
                       colours=c("darkred", "white"), guide="colourbar") + 
  #scale_fill_gradient2(name="LISA p", low="blue", mid="white", high="red", midpoint=0) + #scale_fill_gradient2 defaults to midpoint=0
  #scale_fill_gradient(name="LISA p", low="white", high="darkred", guide="colourbar") + 
  labs(x="", y="", title="LISA p-values (residuals)") +
  theme(axis.ticks.y = element_blank(),axis.text.y = element_blank(), # get rid of x ticks/text
        axis.ticks.x = element_blank(),axis.text.x = element_blank(), # get rid of y ticks/text
        axis.line.x = element_blank(), axis.line.y = element_blank(), # get rid of axis lines, if present
        plot.title = element_text(size = 10, lineheight=.8, face="bold", vjust=1), # make title bold and add space
        legend.title = element_text(size = 8), # legend title size
        legend.text = element_text(size = 8), # legend label size
        panel.grid.major = element_blank(), panel.grid.minor = element_blank(), # removes background grid
        panel.background = element_blank(), # removes grey background; could add black axis lines with #axis.line = element_line(colour = "black")) 
        panel.spacing=unit(c(0,0,0,0), "lines"),
        plot.margin=unit(c(0,0,0,0), "mm")) + # sets margin around full plot at top, right, bottom, and left; units can also be "lines" or "cm"
  coord_equal(ratio=1) +
  geom_polygon(data=shp.st, aes(long,lat, group=group), colour="grey50", fill=NA) 
plot5pyperma

#tiff(file="./figures/figA_lisap-pyperm-res-color.tif", height=6, width=9, units="in", res=300) # could also do pdf, jpeg, bmp, tiff
#print(plot5pyperma)
#dev.off()

# LISA clusters - residuals using permutation method - color

plot6pyperm <- ggplot(shp.df) + 
  aes(long,lat,group=group,fill=as.factor(ls.py.qresp)) + 
  geom_polygon(colour="transparent", fill="white")
#plot6pyperm

plot6pyperma <- plot6pyperm + 
  geom_polygon(aes(x = long, y = lat, group = group), data = shp.df) +            
  #             colour = 'white', fill = 'black', alpha = .4, size = .3) +
  #scale_fill_distiller(name=NULL, palette = "Greys", trans="reverse", breaks = pretty_breaks(n = 5), na.value="white", guide="colourbar") + 
  scale_fill_manual(name="LISA cluster", values=c("white", "red", "blue","lightblue","pink"), breaks = c("0", "1", "2", "3", "4"), 
                    labels=c("n.s.", "high-high", "low-low","low-high","high-low"), guide="legend") + 
  labs(x="", y="", title="LISA Clusters (residuals)") +
  theme(axis.ticks.y = element_blank(),axis.text.y = element_blank(), # get rid of x ticks/text
        axis.ticks.x = element_blank(),axis.text.x = element_blank(), # get rid of y ticks/text
        axis.line.x = element_blank(), axis.line.y = element_blank(), # get rid of axis lines, if present
        plot.title = element_text(size = 10, lineheight=.8, face="bold", vjust=1), # make title bold and add space
        legend.title = element_text(size = 8), # legend title size
        legend.text = element_text(size = 8), # legend label size
        panel.grid.major = element_blank(), panel.grid.minor = element_blank(), # removes background grid
        panel.background = element_blank(), # removes grey background; could add black axis lines with #axis.line = element_line(colour = "black")) 
        panel.spacing=unit(c(0,0,0,0), "lines"),
        plot.margin=unit(c(0,0,0,0), "mm")) + # sets margin around full plot at top, right, bottom, and left; units can also be "lines" or "cm"
  coord_equal(ratio=1) +
  geom_polygon(data=shp.st, aes(long,lat, group=group), colour="grey50", fill=NA) 
plot6pyperma

#tiff(file="./figures/figA_lisacl-pyperm-res-color.tif", height=6, width=9, units="in", res=300) # could also do pdf, jpeg, bmp, tiff
#print(plot6pyperma)
#dev.off()


#################################################
#
# combined graphs Figure 3 (res) in color 
#
#################################################

tiff(file="./figures/appendix/fig3A_comb_res_col.tif", height=9, width=6, units="in", res=300) # could also do pdf, jpeg, bmp, tiff
grid.arrange(plot4a, plot4bpyperma, plot5pyperma, plot6pyperma, ncol=1)
dev.off()

#end